Collisionless shock formation, spontaneous electromagnetic fluctuations 

and streaming instabilities 
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Collisionless shocks are ubiquitous in astrophysics and in the lab. Recent numerical simulations 
and experiments have shown how they can arise from the encounter of two collisionless plasma 
shells. When the shells interpenetrate, the overlapping region turns unstable, triggering the shock 
formation. As a first step towards a microscopic understanding of the process, we analyze here in 
detail the initial instability phase. On the one hand, 2D relativistic PIC simulations are performed 
where two symmetric initially cold pair plasmas collide. On the other hand, the instabilities at work 
are analyzed, as well as the field at saturation and the seed field which gets amplified. For mildly 
relativistic motions and onward, Weibel modes govern the linear phase. We derive an expression 
for the duration of the linear phase in good agreement with the simulations. This saturation time 
constitutes indeed a lower-bound for the shock formation time. 



I. INTRODUCTION 

Colliding plasma shells are present in a variety of phys- 
ical settings. Astrophysical jets produced by black holes 
are expected to generate a shock when interacting with 
the interstellar medium [H, [Ij. Still in astrophysics, the 
Fireball scenario for Gamma- Rays-Bursts [1, 0] relies on 
shock particle acceleration 0-0|: where the shock arises 
from the encounter of two ultra-relativistic plasma blobs 
ejected from a central engine. Non-rclativistic Super- 
nova Remnant Shocks are also instrumental in accelerat- 
ing high energy cosmic rays 0, ^ . 

For the collisionless environments considered, Particle- 
In-Cell (PIC) simulations are an efficient tool to study 
these processes. The formation of a shock following the 
collision of two plasmas was first explored in Ref. 
and observed in Ref. [llj. Subsequent particle accelera- 
tion has been observed in numerous simulations 
In addition, shock generation through counter-streaming 
plasmas has already been observed in laboratory [Tsl - 
|22|, and the conditions required to drive near relativis- 
tic shocks have been identified [1^. The corresponding 
particle acceleration could be a promising alternative to 
current plasma accelerator schemes [24| . 

Although the full shock formation process has thus 
been now repeatedly observed, a first principle under- 
standing of the very birth of the shock is still lacking. 
Such an understanding could provide an accurate timing 
of the shock formation time, and constraints the condi- 
tions required to form a shock in the first place. Whether 
they are in the lab, in a computer or in the vicinity of 
a supernova, it should be possible to separate the sce- 
nario leading to the shock into two phases represented 
schematically on Figure [TJ In the first phase, plasma 
shells make contact, then overlap, and the overlapping 
region turns unstable. An instability grows and satu- 
rates. At this junction, the total density in the overlap- 
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FIG. I: Phases of the shock formation. Two identical pair 
plasmas interpenetrate. The overlapping region turns unsta- 
ble, and two shocks form near the border of each shell. The 
simulation box contains half of the system. 



ping region is roughly the sum of each plasma density. A 
second phase is therefore needed during which nonlinear 
processes pick-up the system from the end of the linear 
phase, and build-up the Rankine-Hugoniot expected den- 
sity jump near the borders of the interpenetrating shells. 

The present paper is concerned with the first of these 
two phases. The collision of two identical cold relativis- 
tic pair plasmas has been simulated in 2D with the PIC 
Code OSIRIS [H, llll . The details of the simulations are 
given in Section [V] This setup has been chosen for its 
simplicity, allowing for a direct comparison with theory 
as the only free parameter is the initial Lorentz factor of 
the shells 70. In the simulation, a neutral e~/e~'" plasma 
is made to bounce back against a wall and to interact 
with itself (Fig. [T]), which enables to describe only half of 
the symmetric physical system. Periodic boundary con- 
ditions are applied in the transverse direction. A series 
of snapshots from a simulation with 70 = 25 is displayed 
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FIG. 2: Integrated density in the direction normal to the flow for 3 instants of a typical shock formation simulation. The last 
plot shows the growth of the magnetic energy integrated over the transverse direction, and xi £ [0, 7y/joc/'^p]- The dashed line 
is the theoretical growth-rate. The initial Lorentz factor was 70 — 25. All the field growth plots look qualitatively the same 
until 70 = 10*. The saturation time Ts is t2. The field at saturation is B{ts). 
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FIG. 3: Same as Fig. El but using 800 particles per cell. 



on Fig. [2] Only the right part of the system pictured 
on Fig. [1] is showed. These successive plots of the in- 
tegrated density along the direction normal to the flow 
clearly show how the overlapping region at twice the up- 
stream density turns unstable, before the shock density 
jump builds up. We observe in the simulations that in 



phase 1 the fields grow in a well defined spatial region 
that extends up to ~ 1 ^f^cjuj^ from the wall. The sat- 
uration time Ts {ti on Figure [2]) is defined as the end of 
the exponential growth of the field energy integrated over 
the region x\ € [0, T^/Toc/wp], where the non-relativistic 
plasma frequency reads uip = Aimee^ /me- The field at 
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saturation is simply B{ts). 

As discussed in Sec. |Vl simulations have been run with 
8 particles per cell. Figure [3] displays the same data as 
Fig. [21 but running the simulations with 800 particles per 
cell. No noticeable qualitative difference appears with 
respect to our runs. 

As will be checked, the instabilities at play can be in- 
terpreted in terms of the homogeneous theory for such, 
although the geometry is finite here, since our shells have 
one contact open boundary. As it amplifies a seed field 
from its initial fluctuation value to saturation, the in- 
stability governs this first phase of the shock formation 
process for a time Ts that we labeled "saturation time" . 
A theoretical determination of this time, in good agree- 
ment with the simulations, is the main result of this pa- 
per. Not only docs give a lower bound for the shock 
formation time, it also sheds a light on the amplitude of 
the amplified initial fluctuations. 



FIG. 4: (Color Online) Growth-rate in ujp units, in terms of 
Z = 'k.vo/ujp for 70 = 1.1 {left) and 70 = 10 {right). 



beams, counter streaming pair beams are linearly equiv- 
alent because the linear regime scales like the square of 
the charge. 

We thus find that unless 70 < ^3/2, current filamen- 
tation should govern the interaction with a growth rate, 



II. INSTABILITY ANALYSIS 

Here we deal with the first phase of the shock forma- 
tion, namely the instability of the overlapping region, 
where we focus on relativistic shocks. Indeed, if counter 
streaming collisionless plasmas were not unstable, they 
would simply go through each other. Let us start ignoring 
the finite geometry at stake here, and analyze the system 
as if it were homogeneous. The full unstable k spectrum 
has been analyzed long ago in the cold regime, where a 
shell is much denser than the other |27h29| . These early 
results were recently generalized to the hot symmetric 
case [3^ [m . For wave- vectors aligned with the fiow, we 
find two-stream unstable modes. For wave-vectors nor- 
mal to the flow, we find the current filamentation, or 
Weibel, instability. Finally, modes propagating at arbi- 
trary angle with the fiow are also unstable. As the two 
plasmas penetrate each other, all the modes are excited. 
But the fastest growing one quickly overcomes the others, 
and shapes the linear phase. For the case we consider, 
a calculation of the growth-rate for every possible wave 
number is pictured on Fig. |4] for two Lorentz factors 
70 = 1.1 and 10; in terms of the reduced wave-vector. 
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where Vq is the initial velocity of the plasmas, and LOp the 
electronic plasma frequency of one of them. The calcu- 
lation, like the simulation, is conducted in the center of 
mass reference frame, where the two plasmas come from 
opposite directions at the same speed. For 70 = 1.1, 
the dominant mode is oblique while for 70 = 10, current 
filamentation dominates. An in-depth study of the prob- 
lem found indeed that only these two types of modes 
can dominate j30l [sij . The transition from oblique to 
filamentation occurs for 70 ~ ■\/3/2, as explained in 
App endix Rl Note that although the analysis of Rcfs. 
[30I Isij was conducted for counter streaming electron 



Comparing this value to the growth of the field observed 
in the overlapping region results in a very satisfactory 
agreement, as evidenced on Fig. [2l We have also checked 
that the Weibel/oblique transition does occur around 
70 — \fij^- Note that a shock also forms for 70 < ■\/3/2 
(not shown). 

This oblique/filamentation transition may seem at 
odds with previous works on instability hierarchy (30l - [3^ , 
suggesting filamentation always governs the spectrum for 
symmetric systems. Electrostatic instabilities with par- 
allel wave vectors have equally been found slower than 
filamentation for relativistic flows [l^H^l- However, the 
relevant hierarchy maps, like Fig. 5 of Ref. [s^ for exam- 
ple, already showed filamentation does not govern sym- 
metric systems all the way down to 70 = 1. Instead, a 
very little gap was found for oblique electrostatic modes 
to dominate, between 70 = 1 and a unspecified value of 
70 > 1. Until now, this little gap did not attract much 
interest, and it is still overall fair to say that in the rela- 
tivistic regime, filamentation is the important instability 
for symmetric systems. 

How can a theory developed for an homogeneous sys- 
tem, apply to the present inhomogeneous system? The 
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instability time scale varies like 8 ( 
time y^^/ujp after contact, the overlapping region is al- 
ready d ^ .^Toc/wp wide. But the parallel scale length 
relative to instabilities is precisely Xi = y/^c/uip. Even 
if at the very beginning of the instability process, d ^ Xi 
is not fulfilled, the strong inequality is quickly realized 
with time passing, so that most of the instability process 
develops in a setting fulfilling the homogeneous approxi- 
mation. 

Knowing the growth-rate ([2]) should allow for an accu- 
rate timing of the linear phase. Assuming the instability 
amplifies a seed field of amplitude Bi up to a saturation 
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level Bs , we can write for the saturation time r. 



Bf^B:e'^' ^ r. 
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where, for convenience, we consider the field energy 
ratio, instead of the field itself. Determining the satm^a- 
tion time amounts then to determine the initial and final 
fields. We will first discuss the saturation field. 



III. FIELD AT SATURATION 




One way to derive the value of the saturation field 
Bf = B{ts), consists in stating that the field grows ex- 
ponentially as long as it is small enough for the system 
to fit the linear approximation. Since a field Bf affects 
particles on a time scale given by the cyclotron frequency, 
this implies [H, HI] , 
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= S ^ 



Bfi = oc. 



(4) 



Another way of evaluating the field at saturation is to 
write that as it grows, particles start oscillating in the 
field of the fastest growing fc,„ at frequency 13, 3^ , 
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Here again, linear theory yielding to an exponential 
growth breaks down when ~ (5, which gives a second 
value for the field at saturation, 
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Finally, one can write the linear approximation breaks 
down when the Larmor radius of the particles in the 
growing field equates This third criterion thus gives, 
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where Vo c has been used. The quantity km has been 
numerically measured with. 
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Accounting in addition for the growth-rate expression 1^ 
gives. 



where, 
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FIG. 5: Field at saturation from the simulations (circles), 
compared with Eqs. (fQl Ulf) . 



is the magnetic field unit of the simulations. 

The magnetic field for the filamentation instability 
grows like sin(fcy)e**. As a result, particles in the vicinity 
of y = [tt] are the ones involved in the second satura- 
tion mechanism, described by Eq. Particles near 
y = TT [tt] experience the kind of trapping involved with 
Eqs. (|4l7p . The linear hypothesis is first broken when the 
field reaches mm{B fi, B f2, B f^) = Bj^. Figure [5] com- 
pares the field observed in the simulation at the end of 
the linear phase with Eqs. (|9llTT|). The agreement with 
Eq. pT|) is good and the correct scaling is recovered. 
At any rate, a numerical pre-factor cannot play a major 
role once inserted into the logarithm of Eq. ^ for the 
saturation time. 

A consequence of the observed 70 scaling is that the 
field energy relative to the beam one reads. 
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displaying the near-equipartition already noted by vari- 
ous authors [l^, HI] . 



IV. THE INITIAL FIELD AMPLITUDE 

We now turn to the evaluation of the initial field ampli- 
tude. The idea is that the instability mechanism picks up 
a seed field from the spontaneous fluctuations of the sys- 
tem, and am plif ies it. Starting with Buneman, Salpeter 
and Sitenko [39l - l41| . various authors have been dealing 
with plasma fluctuations (42| - |48j . 

Though initially focused on stable systems, fluctuation 
theory has been recently extended to weakly amplified 
modes [46l - l48j . Note also that the ability of PIC simula- 
tions to correctly render them has been checked [i^ . 

A first question to ask could be the following: should 
we consider the instability starts from the fiuctuations of 
one single plasma, or from the fluctuations of the system 
formed by the two overlapping plasmas? In other words, 
should we consider the fluctuations of the system before 
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it turns unstable, or after? We will now argue that we 
consider the fluctuations of the stable, isolated plasma 
shells, before they interpenetrate. 

Before they overlap, each plasma shell comes with its 
own fluctuations. As they start to overlap, the fluctua- 
tion fields for each plasma will adapt to each other. But 
on the very same time scale, the instability process be- 
gins. We thus consider the seeds for the instability are 
the ones which were already present in the system before 
the plasmas started to overlap. Filamentation for ex- 
ample, needs unbalanced counter-streaming currents to 
start growing. As they approach each other, both plasma 
shells already display spontaneous fluctuations normal to 
the drift. When they start to interpenetrate, these fluc- 
tuations instantaneously provide the needed unbalanced 
currents to destabilize the whole system. Hence, their 
amplitude will be the amplitude they had before they go 
unstable. 



where the sum runs along the 2 species involved, namely 
electrons and positrons. Both species are assumed to 
obey Maxwell- Jiittner distribution functions [HI, , 
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where K2 is the Bessel function of the second kind and 
vqs the drift velocity of the s-th species. The evalua- 
tion of these quadratures by means of complex analysis 
techniques is extensive and quite involved, and will be 
reported in details in a separate paper [sl] . 

The w-integrated fluctuation energy density reads. 



(20) 



A. Fluctuation power spectrum 

We are interested in the magnetic fluctuation spectra 
of a relativistically drifting, stable plasma. For the non- 
rclativistic case, such calculation has been performed by 
Yoon [11]. For the relativistic case, the amplitude of 
the spontaneous magnetic fluctuations for fc|| = can 
be deduced from the linearized Vlasov-Maxwell system 
which yields the relation between the electric and current 
density fields [50| . 



j(k,w) = Z-i(k,w) • E(k,w), 
where is the tensor 
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where it is assumed that the plasma drifts along the x- 
axis and the wave number fcj^ is along the y-axis. There 
follows [ill, [13, 



(16) 



where EE^ ^.^ is the fourier transform of the sponta- 
neously emitted electric field fluctuation tensor and f 
designates the hermitian adjoint. Taking the xx com- 
ponent of Eq. (|16p and dividing it by the square of the 
phase speed gives: 
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The dielectric tensor eap is given by [5; 

£0/3 = Sap 



(18) 
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Since filamentation modes have w = 0, it is interesting to 
consider the spectrum density for a; = 0. In the regime 
1 <C 7n <g; we consider here, and for k± = Up/c^/jQ, it 
reads 1531, 
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Equation pO|) needs integration over a k domain to ob- 
tain the available power in the corresponding fluctua- 
tions, while result (j2ip needs integration over a given 
(k, uj) domain. 

One could assume the instability process is indeed a 
initial condition problem, so that it can discriminate the 
unstable k's, but not their frequency. In such case, Eq. 
(|20| should be used to derive the initial fleld amplitude. 
We here argue instead that the instability process can 
discriminate the fluctuation frequency, but up to a pre- 
cision ±(5. In other words, cj = is selected for growth, 
but this selection should be inaccurate to an order ±6 be- 
cause during the first growth period, the plasma cannot 
discriminate waves varying at w = ± (5 from the ones 
at = 0. The two approaches will be later compared, 
and the "a; = selection" will be found in slightly better 
agreement with the simulations. 



B. The k-integration domain 



Whether we use Eqs. (|20p or (f2T|) for the available 
energy for growth, we thus need to integrate over a k 
domain. As indicated on Fig. |4l wave-vectors selected 
for growth form a narrow band around the normal axis, 
extending to infinity and of width 



2 



p 



70 c 



(22) 



6 



in the parallel direction (we set here vq ^ c). 

Regarding the integration domain in the normal direc- 
tion, it has already been mentioned that the fastest grow- 
ing mode has been numerically found for k± ^ Wp/c^To. 
We shall thus integrate Eqs. (|20l2ip from k±^min to 

k±,max with, 



1 UJr, 
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where the factors 1/2 and 2 have been arbitrarily chosen 
to bracket k„i = ojp/c^/^. Note that the end result 
is almost independent of these constants because of the 
logarithm function in Eq. ([3]). 



C. Saturation time from oj-integrated fluctuations 

The w-integrated energy density (|20p is now integrated 
in the following way. 
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Clearly, it is an averaged initial amplitude over the modes 
likely to grow the most. The result reads. 
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Inserting this result in Eq. ^ for the saturation time, 
we find 
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, for 70 < ^, 



where n is the plasma density. 



D. Saturation time from fluctuations near uj — 

While the k-integration domain remains unchanged, 
an w-integration domain is now needed. As stated ear- 
lier, physical reasoning would suggest an integration over 
[—(5,(5], because the instability mechanism cannot dis- 
criminate fluctuations with lo = from the ones with 
—S < UJ < 6 during the first growth period or so. Once 
a given fluctuations has been significantly amplified, i.e, 
has grown during ^ (5~^, it will keep on growing. 

But the integration domain is eventually much smaller 
than [—(5,(5] because the density Bk^.u is extremely 
peaked around uj = with. 
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At k±_ = Wp/cyTo, the peak width Suj is given by [Si 

c UJr, 
OLU 



(28) 



which turns out thinner than the growth rate 5 = 
Wp-\/2/7o, especially for /i ^ 1. The energy density (PT|) 
is therefore integrated over [—6uj, 6uj], and the initial field 
amplitude Bi from the energy available for growth com- 
puted as. 



Bl 
Stt 



27rfc I dk 



dku 



BL.o 



5u 

dw 

-&u> 87r 



(29) 

Given the narrowness of the (k, uj) integration domain, 
we set Bk^^ ^ ^fc^.o as given by Eq. (PT|) in the inte- 
grand. A little algebra gives 
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and the saturation time can be cast under the form 
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Finally, it is to be reminded that our theory has been 
implemented for a 3D geometry, whereas simulations are 
2D. The corresponding 2D saturation time is derived in 
Appendix IB] and reads. 
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where N is the number of macro-particles per cell. 



V. COMPARISON WITH SIMULATIONS 

In order to test the theory in the early stage of 
shock formation, ab-initio Particle-In-Cell simulations 
have been performed [l^ll^. The shock is launched with 
the piston-wall method, where two counter-propagating 
symmetric plasma beams are produced by injecting one 
beam by a cathode from one side of the simulation 
box, and being reflected at the opposite wall. Here 
we simulate pair plasmas with low temperature param- 
eter ^ = mc^/kBT = 10^70, and Lorentz factor 70 € 
[25,10^], so that 1 < 70 < /i is fulfilled. The parti- 
cles are injected along the x axis with a temporal res- 
olution At = 0.025y/^/ujp and the size of a cell being 
Ax ~ O.OS-y/Toc/wp, using quadratic interpolation and 8 
particles per cell. Note that this results in a grid size 
much larger than the Debye length, which could trigger 
the grid instability [H^j- However, for the combination of 
grid sizes and number of particles per cell we have used, 
this instability has a much longer time scale than the 
typical times analyzed in our simulations. Moreover, in 
our simulations we use higher order particle shapes and 
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current smoothing (a 5 pass binomial smooth is used). 
This improves significantly the energy conservation prop- 
erties of the algorithm and slows down even further nu- 
merical heating [2^. The present spatial resolution has 
been tested in many shock simulations confirming its ap- 
plicability. Convergence tests have been performed with 
smaller cell sizes, and no significant deviations was found. 
Physically, any numerical heating arising from the grid 
instability is irrelevant to the dynamic of the problem be- 
cause the free energy in the system converted to thermal 
energy is much larger than the energy associated with 
numerical heating. 

We ran test simulations with up to 800 particles per 
cell and observed only small deviations to the results 
reported here (see Fig. [3]). Detailed convergence tests 
were also performed. The two-dimensional box with 
Lx ~ 125.y7oc/ajp and Ly = b^/^c/up has absorbing 
boundaries for the particles along x and is periodic along 
y. For the fields, conducting boundaries are used at the 
perfectly reflecting wall and open boundary conditions at 
the cathode. 

Since we are interested in the early stage of shock for- 
mation, the question arises as to whether the piston- wall 
method is appropriate. We first simulate the periodic 
system of counter-streaming beams, corresponding to the 
model of unstable fiuctuations that is the basis of our 
theoretical approach. In this case, no shock is formed 
and we are able to identify the growth rate and satura- 
tion time of the magnetic field energy. We compare the 
periodic system with the piston-wall setup and, further- 
more, with the full shock formation process, where in x 
direction absorbing boundaries have been used for the 
particles and conducting boundaries for the fields. In the 
latter case, two symmetric shocks are propagating out- 
wards and this allows us to identify non-physical fields at 
the reflecting wall in the piston-wall setup. 

Fig. |6] shows the evolution of the magnetic field energy 
tB normalized by the kinetic energy in the box at time 
zero eo for the three different setups for 70 = 25. The 
comparison shows that the growth rate and saturation 
level of the field is independent of the setup. The theory 
of the periodic system applies to the non-periodic system 
as well, where the overlapping beams go unstable, and 
the fields at the refiecting wall do not seem to affect this 
process. There is only a small deviation in the initial 
fiuctuation level, which for 70 = 25 leads to a shift of the 
saturation time ~ between the different setups. On 
the interesting time scales for the saturation time (see 
Fig. [7]) this deviation is negligible, so that we conduct 
the simulations with the piston-wall setup in order to 
save simulation time. 

Theoretical results are now bridged setting 

\^pJ 0.053 70^^ 

Figure [7] compares the saturation time measured in the 
simulations with Eq. pil) accounting for fiuctuations 
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FIG. 6: Magnetic field energy evolution for different simu- 
lation setups and 70 = 25. Black: piston-wall method, Red: 
full shock picture. Blue: periodic system of counter-streaming 
beams. Black dashed: theoretical growth rate. A detailed de- 
scription of the models is given in the text. 
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FIG. 7: Saturation time TsUJp from the PIC simulations, cir- 
cles, from the fiuctuations near u = Q Eq. (|31[) . bold line, from 
the cj-integrated fluctuations Eq. (|26l) . thin line, and from the 
2D formula (|32p . thin dashed line. The 3D a;- integrated and 
the 2D theories give almost the same result. 



near w = 0, Eq. ([26)) accounting for w-integrated fluctua- 
tions and the 2D formula ([5^ . As expected, considering 
only the fiuctuations around a; = yields a larger satu- 
ration time, arising from a lower initial noise amplitude. 
The slight underestimation of the simulation results can 
be attributed to at least two factors. On the one hand, 
Eq. (|26p necessarily remains a lower limit, as the integra- 
tion domain only brackets the mode selected for growth. 
On the other hand, it is difficult to model the level of fluc- 
tuations in the simulations realistically, since it is depen- 
dent and sensitive on the choice of numerical parameters 
of the simulations. 
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VI. CONCLUSION 



Appendix A: Oblique to filamentation transition 



In this paper, we have analyzed in detail the first part 
of a coUisionless shock formation process. Wc chose the 
simplest possible setup for analytical calculations. Wc 
have thus run relativistic PIC simulations of two inter- 
penetrating cold pair plasmas shells, the initial Lorentz 
factor being the only varying parameter. In the present 
coUisionless conditions, these shells could simply pass 
through each other. But because the overlapping region 
is unstable, turbulence is triggered which eventually leads 
to the shock formation. 



It can be seen from Fig. IHthat the growth rate at large 
reaches a limit bz^.oa which is function of Z\^. For 
7o — 1.1, (5zj^.oo(^||) reaches an extremum for Z\\ ^ 0, 
which corresponds to a spectrum governed by oblique 
modes. Then, for 70 = 10, the extremum is reached at 
Z|| = 0, and filamentation dominates. The first derivative 
dSz^.oo/ dZ\\ always vanishes for Z|| = 0. The transition 
from one regime to the other occurs then when the second 
derivative vanishes at Z|| =0. 

The asymptotic dispersion equation for Z±_ = 00 can 
be determined and reads, 



The shock formation time has been determined from 
the expression of the dominant growth rate, the field at 
saturation and the seed field amplified by the instability. 
The dominant growth-rate could be determined from the 
theory derived for an homogeneous, infinite system, in 
spite of the limited extension of the overlapping region. 
The agreement between the simulations and our simple 
model is due to the fact that although finite, the center of 
the unstable region, where modes start growing, satisfies 
the homogeneity criterion. 

The field at saturation is correctly given by any of the 
3 existing saturation criteria, as the wave number of the 
dominant unstable mode precisely adapts for these crite- 
ria to converge (up to a numerical constant). The reason 
why the system "chooses" to amplify preferably this k± , 
in spite of the absence of a peak in the growth-rate curve 
6{k±) could be the topic of future works. 



4(1 - 7o') - 2{x' + Zjho + [x^ - Z\f^l = 0. (Al) 

This equation can be solved, and the growth rate for 

Z \ = 00 is, 



4, 



y^l + 473(^2 +/3270) 



(A2) 



Deriving twice the expression above with respect to 
Z\\ gives the Lorentz factor for the transition from the 
oblique to the filamentation regime. 



70 



V3 



1.22, 



0.57 . 



Finally, the seed field Bi which the instability picks 
up for amplification is computed from the amplitude of 
the spontaneous fluctuations of one single relativistically 
streaming shell. On the one hand, this density is inte- 
grated over the k domain likely to grow. On the other 
hand, we have tested the Bi value obtained assuming 
the instability mechanism purely acts as an initial value 
process, unable therefore to discriminate the frequency 
of the noise, and the Bi value obtained assuming the in- 
stability selects for amplification those fluctuations with 
cj = 0. By computing the saturation time given by both 
options, we find the second one fits slightly better the 
simulations. 



Appendix B: Application to a 2D PIC plasma 



Care must be taken when using the formula ([2T|) for a 
2D PIC-modeled plasma. The plasma is then composed 
of macro-particles with charge and mass equal, respec- 
tively, to = W^q and Mp = Wpm, where q and m 
denote the real particles' charge and mass, and Wp is the 
statistical weight. In a 3D plasma, Wp is a dimension- 
less quantity, whereas it is a lineic density in 2D. For the 
numerical plasma to behave collectively as its physical 
counterpart, the plasma frequencies of the two systems 
must be equal, which implies 



The reasoning used to time this first phase can in prin- 
ciple be adapted to any settings. By the end of the linear 
growth phase, the density of the overlapping region is still 
about twice the upstream density, as evidenced on Fig. 
[2I Indeed, because the linear regime requires small per- 
turbations, it is necessarily over by the time density per- 
turbations reach Sn ^ 0.1 — 1. For the present system, 
the density jump around the shock soon to be formed 
is around 3. This implies other processes have to pick- 
up the system from the saturation time up to the shock 
formation. We plan to dedicate future works to these 
mechanisms. 



mujp AxAy 
47rg2 N 



(BI) 



where N is the number of macro-particles per cell and 
Ax ~ Ay is the cell size. 

In a 2D geometry, the fluctuation field is then given by 



Bf 



\,min J —k\\ 



12 70 2 

— WpTnc 
n \ c J 



dku 



r2 



(B2) 
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Note that the normalized inverse temperature /i is also 
an invariant. Substitution of Eq. (jBl[) and Aa; = 
0. 05 y^c/cjp into (jB2p readily yields 
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2.5 X 10^3 |3 ^2 /^cwp 



TV 



There follows the ratio 



S,^ V 3 70 



and the saturation time given by Eq. p2p . 
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